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We present a Hamiltonian in real space which is well suited to study numerically the behavior of 
holes introduced in III-V semiconductors by Mn doping when the III 3+ ion is replaced by Mn 2+ . 
We consider the actual lattice with the diamond structure. Since the focus is on light doping by 
acceptors, a bonding combination of III and V p-orbitals is considered since the top of the valence 
band, located at the T point, has p character in these materials. As a result, an effective model 
in which the holes hop between the sites of an fee lattice is obtained. We show that around the 
F point in momentum space the Hamiltonian for the undoped case is identical to the one for the 
Luttinger-Kohn model. The spin-orbit interaction is included as well as the on-site interaction 
between the spin of the magnetic impurity and the spin of the doped holes. The effect of Coulomb 
interactions between Mn 2+ and holes, as well as Mn 3+ doping are discussed. Through large-scale 
Monte Carlo simulations on a Cray XT3 supercomputer, we show that this model reproduces many 
experimental results for Gaj -^Mn^As and Gaj-^Mn^Sb, and that the Curie temperature does not 
increase monotonically with x. The cases of Mn doped GaP and GaN, in which Mn 3+ is believed 
to play a role, are also studied. 

PACS numbers: 71.10.Fd, 74.25. Kc, 74.81.-g 



I. INTRODUCTION 

The discovery of high Curie temperatures (Tq) m some 
diluted magnetic semiconductors (DMS) 1 -' 2 -' 3 - has created 
a renewed interest in the study of the properties of these 
materials due to their potential role in spintronics devices 
if above room temperature Tc's are achieved^ 

When a III-V compound is magnetically doped the 
added holes can gain kinetic energy by moving in the 
bands formed by the hybridized orbitals or magnetic en- 
ergy by interacting with the spin of the localized mag- 
netic impurities. The kinetic term induces derealization 
and it is usually studied in momentum space; the mag- 
netic term tends to localize the carriers and it is natural 
to study it in real space. In the limit of large Mn doping 
the kinetic term should prevail while in the very diluted 
regime the magnetic term has to be accurately consid- 
ered. Most studies of DMS are based in one of these 
two extreme limits. However, in order to understand 
the intermediate region between both extremes powerful 
non-perturbative techniques are needed. The goal of this 
work is to provide such a tool. 

The band structure of Zinc-Blende type semiconduc- 
tors has been accurately obtained using a variety of differ- 
ent approaches such as pseudopotential methods, tight- 
binding techniques, and by careful consideration of the 
symmetries involved through the application of group 
theory. 5 . However, many of these calculations involve a 



large number of orbitals per site making them unsuit- 
able for simulations of doped systems with present day 
computers. The successful Luttinger-Kohn model^ that 
describes the top of the valence band of semiconductors 
with diamond structure is very useful to study hole dop- 
ing in cases of strong hybridization. When the doped 
holes are not localized they go into the valence band and 
a formulation of the problem in momentum space is ap- 
propriate. This approach, known as the valence band or 
hole-fluid scenario, has already been applied to DMSZ^ 
providing results that would be valid in a regime in which 
the doped holes are uniformly distributed in the crystal 
and the randomness in the impurity distribution becomes 
irrelevant. On the other hand, it is well known that in the 
limit of strong dilution the doped holes will have a ten- 
dency to become localized at the Mn impurities. In this 
extreme regime holes move by hopping through overlap- 
ping orbitals and this behavior is better described in real 
space. The study of this limit has been carried out only 
at a phenomcnological level using tight binding Hamil- 
tonians that, while taking into account the random po- 
sition of the impurities, do not provide a realistic band 
description of the bands resulting from the overlaps of 
the atomic orbitals involved ^ 10 ' 11 : 12 ' 13 : 14 

The goal of this paper is to provide a realistic, yet sim- 
ple, tight-binding-likc Hamiltonian, in real space, that 
reproduces the top of the valence band of Zinc-Blende 
type semiconductors with the smallest possible number 
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of degrees of freedom per site, and that handles the in- 
teractions between the doped holes and the randomly 
distributed magnetic impurities. It will be shown that 
six orbitals per site are needed to capture experimental 
properties of the DMS even in the case of a strong spin- 
orbit (SO) interaction in which, naively, only four orbitals 
are relevant at the top of the valence band£ This model 
reproduces the Curie temperatures observed experimen- 
tally for various values of Mn doping and compensations 
and it captures the Curie- Weiss behavior of the magne- 
tization curves in the low compensation regime. It also 
predicts that for Mn doped GaAs the Curie temperature 
would pass through a maximum value Tq ~ 220K for 
x » 12% suggesting that higher Mn doping is not the 
route to achieve room temperature Tc with GaAs. 

The paper is organized as follows: the general tight 
binding approach is described in Section II; Section III 
is devoted to systems in which the spin-orbit interaction 
can be neglected; the spin-orbit interaction is considered 
in Section IV by working in a basis formed by total angu- 
lar momentum eigenstates \ j, rrij >; Section V is devoted 
to the discussion of Coulomb interactions and how to 
handle the possibility of Mn 3+ doping. Numerical results 
in finite systems for a variety of materials are presented 
in various subsections of Section VI; conclusions and a 
summary appear in Section VII. In the Appendix the 
change of basis matrices are provided. 



II. TIGHT BINDING APPROACH 

Although we will provide a Hamiltonian that can be 
applied to any III-V semiconductor with Zinc-Blende 
type structure, we will use Mn doped GaAs as our ex- 
ample. In GaAs, the Ga and As atoms bond covalently. 
Both atoms share the electrons that they have in their 
4s and 4p shells. Ga shares 3 electrons while As shares 
5. The hybridized orbitals have character sp 3 . Although 
the s and p orbitals have to be considered in order to 
obtain the correct band structure, we are interested in 
light hole doping for which we will focus only on the va- 
lence band, since to facilitate numerical simulations we 
need to minimize the number of degrees of freedom. It 
is well known that around the T point the valence band 
of GaAs has j = | character, which arises from the orig- 
inal p orbitals.- Thus, in order to construct the simplest 
model that captures this feature we will consider only the 
p orbitals. The three p orbitals p x , p y and p z in each ion 
can be populated with particles with spin up or down. 

To study GaAs we should consider two interpenetrat- 
ing fee lattices separated by a distance (ao/4, ao/4, ao/4) 
were ao is the underlying cubic lattice parame- 
ter. The Ga ions seat in one of the fee lat- 
tices and the As ions in the other. Each ion 
has 4 nearest neighbors of the opposite species 
located at (a /4, a /4, a /4), (a /4, -a /4, -a /4), 
(-a /4, a /4, -a /4) and (-a /4, -a /4, a /4). 

Since we are only interested in obtaining the valence 



band we are going to consider the bonding combinations 
of the Ga and As p orbitals i 5 ' 15 ! 16 This leads to an effec- 
tive fee lattice with three p bonding orbitals at each site 
that can be occupied by particles with spin up or down. 
Thus, working in this \p, a > basis there are 6 states per 
site of the fee lattice. We will consider the nearest neigh- 
bor hopping of these particles to construct the effective 
tight-binding Hamiltonian. The twelve nearest neigh- 
bors are located at (±ao/2, ±ao/2, 0), (±ao/2, 0, ±ao/2), 
(0, ±a /2, ±ao/2) considering the four sign combinations 
for the three sets of points provided. 

In order to calculate the hoppings we follow Slater and 
Koster.— The nearest neighbors in our effective fee lattice 
are the second nearest neighbors in the original diamond 
structure. From Table I in Ref[T3 we see that the relevant 
overlap integrals in this case are 

E xx = l 2 (ppa) + (l-l 2 ){pp7r), 
E xy = lm[(ppa) - (jppir)], 

and 

E xz = ln[{ppa) - (ppir)}. (1) 

The 12 nearest nei ghb ors in the fee lattice are labeled 
(p, q, r) following Ref 1 1 71. For this geometry, two of the 
indices taking the value ±1 and the remaining one is 0. 
Since 1 = / I V 2 i 2 with b'9' r ) = (M,0), etc. (we 
are following Slater's notation), it follows that I, m, and 
n are equal to or ±l/\/2. Then the hoppings to the 
twelve neighbors that we will label by (/i, v) with and 
v taking the values ±x, ±?/,and ±z, are: 

-C = E ™ (M, v) = 2 [(P?*0 + (PP*)] = -*L (2a) 
if cither /i = a or v = a or 

-C = E xx {^ v) = (ppTc) = -fcL (26) 
if neither v nor \i are equal to a, and 

-*ab = ExyiV, v ) = ±\[(jPPv) - (PP*)} = Tt xy , (2c) 

with the minus (plus) sign for the case in which \i and 
v have the same (opposite) sign. Also notice that the 
interorbital hopping is only possible when (/i, v) and ab 
are in the same plane, i.e., there is no perpendicular in- 
terorbital hopping. 

III. THE MODEL WITHOUT SPIN-ORBIT 
INTERACTION 

If the doped Mn ions go into the parent compound as 
Mn 2+ replacing Ga 3+ , it means that a hole is doped into 
the p hybridized orbitals which form the valence band 
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of the undoped semiconductor. This allows us to write a 
Hamiltonian that takes into account the nearest neighbor 
hopping of the holes in the p orbitals using the hoppings 
calculated above, and their magnetic interaction with the 
localized randomly doped Mn 2+ ions: 

H = \ 12 C ( C U,a C i+P+-«,b + h.C.) 
i,/x,z/,a,a,b 

+J^si a -Si, (3) 

I,a 

where c! creates a hole at site i = (i x , i y , i z ) in orbital 

a with spin projection a, si a =^Q/3 4, a!) a (J «/3 c i./3.a is tlic 
spin of the mobile hole, the Pauli matrices are denoted 
by a, Si is the localized Mn 2+ spin 5/2 at site I (covering 
only a small fraction of the total number of sites N since 
Mn replaces a small number of Ga). t\" b are the hopping 
amplitudes for the holes that were defined in Section II, 
and J > is an antiferromagnetic (AF) coupling between 
the spins of the mobile and localized degrees of freedom. ^ 
The density (n) of itinerant holes is controlled by a chem- 
ical potential. The sites i belong to an fee lattice and the 
versors [iv indicate the 12 nearest neighbors of each site 
i by taking the values ±x, ±y, and ±z, with [i ^ v. 

Notice that we have only three different hoppings: two 
intraorbital ones and via. and the interorbital ones i a t,. 
All have the same absolute value (but not always the 
same sign) for all combination of orbitals and neighbors. 
The interorbital hoppings that have the sign reversed are 



the ones towards sites labeled by [i and v with opposite 
signs, i.e. —t%£ = ,y - Also notice that the interorbital 
hoppings occur only in the planes defined by the two 
orbitals, i.e., they vanish in the direction perpendicular 
to the plane were the two orbitals are = 0. This can 
be seen in the expressions provided in Eq.(l). 

In order to obtain the material specific values of the 
hopping parameters we will write the Hamiltonian matrix 
in momentum space for the undoped case, i.e., Eq.(3) 
with J — 0. We can use Table II or III in RefUJ to do 
this task. The result is given by, 

(Tx txySxy txy^xz \ 

txy^xy Ty txySyz I 3 (4) 

txy^xz txySyz T z I 

for spin up and an identical block for spin down. Here 

T x ^txx(CxCy ~\~ C x Cz) -\~ At xx CyC Z: Ty ^txx{c x Cy ~t~ 

c y c z ) + At xx c x c Zl T z = 4tL {c x c z + c v c z ) + At xx c x c y , and 
Sij = AsiSj with Cj = cos(afci), Sj = sin(afcj), where 
in Slater's notation a = ao/2, with ao the Zinc-Blende 
lattice constant and ki are the momentum components 
along i = x,y, or z. 

We mentioned in the introduction that Luttinger and 
Kohn studied the movement of holes in the valence band 
of semiconductors with the diamond structure. Working 
in momentum space they found an expression for the 
Hamiltonian matrix that describes the top of the valence 
band, i.e., the neighborhood of the T point. In the \p, a > 
basis the matrix has the form: 



Ck x ky 

(2 k'r. 



k 2 z ) 



Ckx ky 
Cky k z 



Ck x k z 
k z ) Ckyk z 

Akl + B{kl + kl) 



(•5) 



There is an identical block for spin down. A, B and C 
are constants that can be defined in terms of the Lut- 
tinger parameters 71, 72, and 73 II3 ) which are material 
specific. The accepted values for GaAs are (71,72,73) = 
(6.85, 2.1, 2.9)££ The constants are given by: 



h 2 

4 = -o~(7i 
2m 



-472), 



Remembering that Eq.(5) is an approximation valid 
for the top of the valence band (i.e. the T point), we can 
expand the cosines and sines in Eq.(4) and we obtain 
the matrix shown in Eq.(5) if we disregard a constant 
term along the diagonal that just shifts the top of the 
valence band to 8tl x + 4$xx instead of 0. Comparing 
the coefficients we obtain expressions for the hoppings in 
terms of the Luttinger parameters and the lattice con- 
stant ao = 2a>^ 



B = - — (71-272), 



^ = -6^73, 
2m 

where m is the mass of the bare electron. 



(6) 



8ma- 



.( 7l + 4 72 ) = 



t 1 - = 



r(7i ~ 872) 



2mag 



2ma\ 



(71 +472), 



(71 - 872), 



4 



_ 3ft 2 
xy Ama 2 



3ft 2 



73 



r73- 



(7) 



Since our goal is to write a tight-binding Hamiltonian 
for holes that will dope at most the bottom of the band, 
we will reverse the signs of the hoppings since the band 
obtained in Eq.(5) gets reflected with respect to zero by 
reversing the signs of A, B, and C— Then, for GaAs, 
ao = 5. 64 A,— and from Eq.(7) we obtain: 



t\ x = -1.82eV, 



1.20eV, t. 



-2.08eV. (8) 



It is important to notice that the evaluation of the 
hopping parameters given by Eq.(7) requires three inde- 
pendent parameters to fix their values. However, if wc 
look at the expression for the hoppings given in Eq.(2) 
in terms of overlap integrals it would appear incorrectly 
as if only two parameters were necessary and the three 
hoppings should be interrelated. The evaluation Eq.(7) 
is more accurate though because it considers, in a phe- 
nomenological way, the influence of the neglected bands 
in the shape of the valence band at T. 



with the diamond symmetry, Luttinger and Kohn showed 
that at the Y point the states with j = | get separated 
from those with j '• = | which are the relevant states at 
the top of the valence band. As a result, naively only 
four states per site, instead of six, become relevant when 
the spin-orbit interaction is considered and the doping 
is light-- However, as it will be shown in Section VI, at 
the levels of doping of interest (x > 3%) the heavy and 
light hole states that become populated are sufficiently 
separated from the Y point that the j — 1/2 contribution 
becomes important. 

To take into account the spin-orbit interaction we will 
have to perform a change of basis from \p, a > to \ j, m,j >. 
This change of basis has been studied by Kohn and 
Luttinger. 6 

The Luttinger-Kohn matrix (Eq.5) in the |j, rrij > ba- 
sis is presented in Eq.(A8) of RefQ. 
The six states that form the basis are: 



1>= 



3 3 
! 2'2 



>, |2>= 



1 



>, 



l»-l|.5>. 



IV. SPIN-ORBIT INTERACTION 

The spin-orbit interaction is important in GaAs and 
in other III-V semiconductors in which the ion V has 
a large mass. It mixes the angular momentum (1 = 1 
for the p orbitals) with the holes' spin degrees of freedom 
producing states with j = § and 3 = \- In a cubic lattice 



M I 3 3 I. I 1 1 I* I 1 1 

|4>=l2>"2>. I 5 >=l2'2 > ' |6>=l 2'-2 > ' 

(9) 

/TO 

Applying the same change of basis to U = ( „ j, 

with T given in Eq.(4) we obtain the 6x6 matrix U' = 
MUM^ 1 where M is the change of basis matrix provided 
in Appendix I. 



/ H hh 


—c 


-b 





b 






—c* 


Hih 





b 


-b\f^ 

V 2 


-d 




-b* 





Hih 


— c 


d 









b* 


-c* 


H hh 




b* 




b* 
V2 


-b,f^ 

V 2 


d* 


-cV2 


H so 







\c*V2 


-d* 


V 2 


b 





H so 


) 



(10) 



with a-a<^2l<\\ ± - t^ x ) 

I i II V ~o~\^x,x — tx,x) c x c y » r~ \ c x c z + CyC z )), 

Hhh = m x c xCy + 2(t£ + tlJ{c x c z + c y c z ), 6 3 V2 



2 " L ^ L 

Hih = ^(S*!,! + t x ,x)( c v c z + c x c z ) + -{t x ,x + 2t XtX )c x c y , 



-4 
V3 



H so = A '— —(c x Cy + c x c z + c y c z ) + A so , (11) 



where Ago is the spin-orbit splitting which is tabulated 
2 || , Ai for different materials^ The matrix in Eq.(10) agrees 

C = 71 ~ x ' x ' {CxCz ~ CyCz > ~ 7/^ x < vSxSv ' with Eq.A8 in Ref 7 in the limit of small ft. 
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A. Hoppings between \ j, irij > States in Real Space To do this, we need to express the c operators that 

appear in Eq.(3) in terms of the c operators in the basis 
Now we need to calculate the hoppings in real space \J, m j > with the help of the basis transformation given 
between the 6 orbitals characterized by j = | and mj = m the Appendix: 
±| and ±| and j = \ with mj = ±|. 



(12) 



and 



M _1 



/ C; 3 3 \ 

*> 2 > 2 

Cj 3 _l 

' 2 ' 2 

C: 3 l 

2 ' 2 

C: 3 _3 
> 2 ' 2 

q I i 

L > 2 ' 2 

Cj i _I 

2 ' 2 



c i,T,z 
c i,J.,ir 
C i>b!/ 



We obtain: 



4 _ 



= cr(—=c: , 3 H =cJ 3 x ) =cl 



and 



V6 



V6 



V3 



/ 2 t At 
W" 1/3^2- 13'*^' 



(14a) 



.-1 1 1 

C( ^C: 3 „3 H = C; 3 1 ) — C: 1 „ 1 . 



cr 5 > CT 2 3 



(15a) 



(13) Now let us consider the inter-orbital hoppings between 

orbitals with j = 3/2: 



,x.z 

era,— era' 



era,— era' 



12 



'•a, -a' — I'a-aV — Vl/2,cr3/2 ~~ \ l <tZ/2,<j1/1> 



V3 ^ 



and 



f x ' z = a v f x ' v 

era, era' ' era,e> 



t M." = n 

''era, — era u i 



0. 



(156) 



where a ^ a' in the above expressions. 

There are no interorbital hoppings between the two 
j = 1/2 orbitals, i.e, 



., 1 1 ia 

Ci.er,v = H^ C i 2 CT 3 H — Cj 3 1 j — C: 1 _ 1 



''era, — era V J 



(15c) 



if a = i, i.e., the symbol that we use to represent mj for 
the orbitals with j = \. 

The interorbital hoppings between j = 3/2 and j — 
1/2 orbitals are given by: 



Replacing in Eq.(3) and rearranging the terms we find 
that the intraorbital hoppings are given by: 



f 11 -4- 9/-L 

,x,y _ ,|| ,x,y _ L xx T ^ xx 



,y,z _ ,x,z _ 

Cr4,cr4 er4,cr4 



(txx + t xx ) 



^1 „ 1 — ''1 ^ 1 

CT 2^2 ^S'^S 



(5t XX + txx) 



A.XZ 

' 1 „3 t _ 3 1 

^2^2 ^^^'^^ 



V6 

a 2>< 7 2 2 ' 2 V6 

1 1 ^ 1 1 y — 

<^2 - CT 2 CT 2 - ff 2 \/2 

V2' 



iff 



er 5 ,-cr 5 ^5 " ff 2 



6 



V2(t± x - tL) 



3\/2 



(15d) 

The actual numerical values are obtained using 
Eq. (7) i 2 — Notice that the hoppings can also be read from 
the matrix elements in Eq.(lO). 



B. Hund Coupling and the Classical Limit of the 
Localized Spin 



The Hund term in Eq.(3) has the form 
Hi = JSi • si, 



(16) 



where Si is a dimensionless quantum spin 5/2 
operator ; 23 ' 2 Si is a dimensionless quantum spin 1/2 op- 
erator, and J, which has units of energy (eV), is deter- 
mined experimentally ^ 18 ' 25 

In order to perform numerical calculations avoiding the 
sign problem 2 ^ which arises when four-fermion interac- 
tions, such as in Eq.(16), are decoupled, we are going to 
take the classical limit for the localized spins. This is a 
good approximation given the large value, 5/2, of this 
spin. The classical limit is given by h — > 0, S — > oo. In 
order to take the limit correctly Eq.(16) is rewritten as 



fli = J[5(5 + 1)] 1/2 5i-si, 



(17) 



where Si is a unit spin quantum operator. 

It is important to realize that the Hund coupling J is 
proportional to h since the spin operators in Eq.(16) are 
dimensionless. This means that lim^o J = 0. However, 
limft_>o,S— ►(» J[S(S +1)} 1 ? 2 = lim R ^o,s-+oo JS, should be 
finite. Defining 



lim JS = J c 

h — >0,S — >oo 



in the classical limit Eq.(17) becomes 
Hi = J c Si ■ si, 

where now 



Si = (sin6*icos0i, sin^isin^i, cos^i), 



(18) 



(19) 



(20) 



(22) 



which is equivalent to 

lim hS = h = 6.58 x lQ- 16 eV, 

i.e., the classical limit according to Refl28l. 

Additional support for these results is obtained mea- 
suring the splitting of the 4 degenerate states at the top 
of the valence band as a function of J. Experimentalists 
use this split to determine Ji 23 i 24 i 29 According to mean- 
field calculations the splitted levels have energy ±Bg and 
±3B G with 8 .* 2 ^ 

Bq = — JxS, 



6 

with S = 5/2. 23 - 24 According to Eq.(23) 



B 



G 



lim Bg = -J c x = - Jx. 
+o,S^oo 6 6 



(23) 



(24) 



Our numerical results indicate that the split of the top 
of the valence band depends linearly on J for values up 
to J < 0.5 eV (see subsection VI. C). In Fig.l we show 
the mean field value of Bg, i.e., Bq F , calculated using 
Eq.(23) as a function of J for x — 0.085 and x = 0.05 indi- 
cated with doted lines. We also display the extrapolated 
linear behavior obtained numerically for J < 0.5 eV and 
the same values of x at T — * with our Hamiltonian. The 
curves indicating Bq in the figure have been obtained by 
monitoring the energy of the four lowest eigenvalues of 
Eq.(31) as a function of J. In the linear regime we ob- 
serve that the ratio between Bq F given by Eq.(23) and 
the measured numerical value in the classical approxima- 
tion Bq is equal to S = 2.5 as expected. The excellent 
numerical agreement between Bq F / S and B C G support- 
ing the result of Eq.(24) is clear in Fig.l. 
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2 3 
J(eV) 



is a classical unit vector, si is still a quantum spin-1/2 
operator and J c is a parameter. The value of J c should 
be determined by comparing numerical results to exper- 
imental data. For Mn doped GaAs comparing numerical 
values of the Curie temperature to experimental data, we 
have observed-! that J c w J. This indicates that 



lim 



JS = J, 



(21) 



FIG. 1: (color online) B^ IF /S vs J for x = 0.085 and 0.05 
obtained using Eq.(23) (doted lines) and Bq vs J calculated 
numerically (extrapolating from the linear regime) with our 
model for the same values of x (dashed lines). 

Now let us consider the Hund coupling term in Eq.(3) 
in the \j, rrij > basis. Eq.(A.lO) in Ref0 has expressions 
for the spin operators in the \j, rrij > basis: 
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(25c) 



Then 



where 



with 



JSi • sj -> J(Sfsf + S*«» + Sf sf ) 



ci 



Hi i 



a t 
,Sj = Ci'SqCi 



ch c\ 



(26) 
(27) 
(28) 



and under the classical approximation that we just dis- 
cussed, 



Sf = sin6>iCOS</>i, S{ = sin^sin^, S( = cos6>;. (29) 

Notice that the Mn ions replace Ga so they will be present 
in a subset of the points of the fee lattice that we are con- 
sidering. The values of J can be obtained from comparing 
experimental data with theoretical models and are ma- 
terial dependent. In the notation of Ref@, J in Eq.(26) 
is given by /3No where /? has units of eVnm 3 and Nq is 
the concentration of cation sites given by 4<2q 3 where do 
is the lattice parameter of the material. (3 is assumed to 
depend only on the characteristics of the parent material 
and it is considered the same for all III-V semiconduc- 
tors. Notice that /3 is called J p d by some authors.- 7 Ref.8 
estimates J for GaN and other materials assuming that 
the accepted value for GaAs is accurate and given by 
J = 4/3<2q 3 = — 1.2eV. x ~ However, in the literature the 



value of J for GaAs ranges between — 0.896^ > J > 
— 3.34eV r > 7 - ' 18 i 25 Using the parameters for GaAs we can 
estimate that — O.OAeVnm 3 > Pm-v > ~ 0.15eVnm 3 . 
Then, J for a general III-V material M is given by 



Jm 



III-V 



(30) 



The sign depends on the definition. It seems to be 
antiferromagnetic for GaAs which means that we need 
to take it as a positive number in our Hamiltonian. The 
sign of J is still under discussion for some of the III-V 
compounds such as GaN^i However, under the classical 
spin approximation discussed above, the sign of J in the 
Hamiltonian becomes irrelevant since the eigenvalues of 
H are the same for both signs. 

Notice that using the parameters provided in this pa- 
per, numerical calculations of the critical temperature 
will be obtained directly in eV. 



C. Six-Orbital Real-Space Hamiltonian 

Combining the results in the previous subsections the 
Hamiltonian in the |j, rrij > basis can be written as: 



H 



o ( t L,Q'b C La Ci +M+","'b + n - C 



+A so V c\ 

— ' l.t 



I 



(31) 



h or |, a and a' can be 1 



where a, b take the values ^ 

or —1, and the Hund term is given by Eq.(26). Notice 
that /j, + v are the twelve vectors indicating the twelve 
nearest neighbor sites of each ion at site i and that I 
are random sites in the fee lattice. To reproduce the 
correct physics we will have to work with an almost filled 
valence band, i.e., a very low density of holes. Eq.(31) 
is the Hamiltonian that has to be programmed in the 
computer simulations. Since there are six states per site, 
a 6N S x 6N S matrix needs to be diagonalized at each step 
of a Monte Carlo (MC) simulation; N s is the number of 
sites in the fee lattice. 

The material dependent hopping parameters are given 
by Eq.(7) and (15), and the Hund coupling by Eq.(30) 
(or from experiments if available). A50 is given by the 
magnitude of the gap between the top of the valence band 
and the split-off band induced by the spin-orbit interac- 
tion. Its value is 0.341eV for GaAs, 0.75eV for GaSb, and 
0.017eV for GaN£ Thus, it is clear that all the parame- 
ters in the proposed Hamiltonian are fixed, i.e., there are 
no free parameters. 
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V. ADDITIONAL EFFECTS OF THE DOPED 
MN IONS 

In this section we will consider, at a phenomenological 
level, the effects of Mn doping that go beyond the mere 
introduction of a localized spin at the doping site I and 
the addition of holes into the system. 

Up to this point we have implicitly assumed that the 
Mn d levels are deep into the valence band of the parent 
compound so that when a Mn ion replaces Ga in GaAs, 
Ga 3+ is replaced by Mn 2+ and a hole is introduced in 
the p orbitals. In this case, negative charge with respect 
to the background at the Mn site creates an attractive 
Coulomb potential for the hole which, in the case of one 
single Mn ion, contributes to producing a bound state at 
O.leV above the top of the valence bandi 24 ' 25 This bound 
state is expected to generate an impurity band at least for 
very light amounts of doping. This impurity band, due 
to hybridization, eventually will merge with the valence 
band as the amount of doping increases. For GaAs it is 
believed that the merging occurs for x < 0.0 1 33 i 34 and 
for this reason the Coulomb interaction is not expected 
to play a role in the relevant doping regime. However, 
the explicit addition of a Coulomb term to the Hamilto- 
nian is important in order to reproduce the GaAs behav- 
ior at very low dopings and, in materials for which the 
hole binding energy is larger, such as GaP and GaN. 35 i 36 
Thus, below, in subsection A, we are going to describe 
how to introduce Coulomb attraction in our model. 

Another effect that is in general ignored in studies of 
(Ga,Mn)As is the presence of Mn 3+ ions. While experi- 
ments appear to indicate that Ga 3+ is replaced by Mn 2+ 
in GaAs ; 23 ' 32 Mn 3+ has been reported in GaP— and 
GaNi 37 ' 38 Thus, in order to extend the present model 
to other materials, in subsection B we will describe a 
phenomenological way of considering the Mn d orbitals. 



A. Coulomb Potential 

The Coulomb potential between the localized impuri- 
ties and the doped holes is neglected in most models for 
DMS^ ' 10 ' 11 ' 12 ' 13 ' 14 since the magnetic interaction appears 
to be sufficient to capture qualitatively many properties 
of these materials including ferromagnetism, and because 
it is believed that at the levels of doping for which the ma- 
terial is metallic it will not play a relevant role^i 3 - 4 How- 
ever, these assumptions need to be confirmed by actual 
calculations. Very few attempts to include the Coulomb 
interaction studying its effects with unbiased approaches 
have been carried out. The case of a single Mn impurity 
in GaAs was studied in Refl39l. They considered the long 
range Coulomb potential suplemented by a central cell 
correction with a square-well shape that is routinely ap- 
plied in calculations of bound state energies for impurities 
in semiconductors^ This is the portion of the potential 
that is expected to be relevant at higher doping since the 
long-range part will be screened. This effect has been in- 



corporated by many author o 3 i i 43 by considering an 
on-site central-cell potential of the form: 

H c = -Vj2 n i- (32) 

The value of V varies widely in the literature. Many au- 
thors determine it by finding the value that added to J re- 
produces the energy of the single hole bound state i 34 ' 39 i 41 
However, the value of V depends strongly on the charac- 
teristic of the model used such as the total bandwidth, 
the number of orbitals, etc. In other investigations V 
is considered a free parameter that can even take nega- 
tive values, i.e., repulsive rather than attractive accord- 
ing to the notation of Eq.(32) i 42 i 43 Repulsive potentials 
are sometimes introduced to reproduce apparent x de- 
pendences of the magnetic coupling J observed in some 
experiments on Mn doped GaN 44 and CdSi 4 ^ 

However, it was pointed out in Refl34 that an on-site 
range potential prevents the overlap of localized hole- 
wave functions that should occur with increasing dop- 
ing. It was proposed that a nearest-neighbor range po- 
tential should be used, and MC calculations performed 
on a highly doped single orbital model showed that an 
impurity band generated by the on-site potential van- 
ished when a nearest-neighbor range potential of the 
same strength was used. In the fee lattice considered 
here, the extended potential has the form 

H 7 6 n = -Vj2(w ni + £ j^—m+v-u), (33) 
1 Nl 7t N w 

where V is added at the site I in which the impurity 
is located and at its 12 nearest neighbors I+/i + v. Nj 
indicates the number of impurity sites that surround the 
impurity site I. Calculations using both Eq.(32) and (33) 
will be presented in Section VI. E. 

B. Mn 3+ 

The phenomenology described in the previous subsec- 
tion is expected to occur when the Mn d orbitals lie deep 
into the valence band so that Mn 2+ replaces the III 3+ 
ion and the doped hole goes into the p orbitals. How- 
ever, if the Mn d levels were in the parent compound's 
gap, the III 3+ ion could be replaced by Mn 3+ with the 
hole in its 3d shell. The 3d shell of Mn is splitted into 
three ti g states and two e g states due to the crystal field. 
It is believed that the hole goes into the tig levels^ 7 - Al- 
though no Mn 3+ has been observed in GaAs, there are 
some indications of it in GaP and GaN^ 7 ^ 

Ideally, at least the three ti g orbitals in each Mn should 
be considered but we would not be able to study the 
problem numerically due to the large number of addi- 
tional degrees of freedom. Thus, we will follow a more 
phenomenological approach by introducing only one ex- 
tra orbital, which will add two extra degrees of freedom 
per Mn due to the spin of the hole. For this purpose, 
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we extend the six orbital model described by Eq.(31) by 
adding an extra orbital at the impurity sites that can be 
populated by holes with spin up or down. To simplify 
the calculations we will consider an extra s-like orbital. 
Thus, a contribution H p _d given by 



(t ss d\ a dj, a + h.c.) + 



X (W4,<r C .7> + h - c -) + J X] 4,7^7,5^,5 ' S I ( 34 ) 



-T,7,<5 



will be added to Eq.(31). In Eq.(34), Vd is the parameter 
that determines the position of the Mn d level in relation 
to the top of the valence band. Its value could, eventu- 
ally, be obtained from experiments but here, in this first 
effort, it is going to be treated as a free parameter; \i is 
the chemical potential that should be added to Eq.(31) to 
fix the number of holes; nj = d\ a dj M is the number 
operator for the holes in the d-orbital; d\ a (dr, CT ) creates 
(destroys) a hole with spin a in the (single) "d" -orbital 
at impurity site /; (/, J) indicate nearest-neighbor impu- 
rity sites; t ss is the phenomenological hopping between 
holes in the d-orbitals of nearest-neighbor Mn impurities 
(notice that this hopping is active only in the disorder 
configurations in which there are Mn ions next to each 
other) . t sp is the hopping between holes in the (i-orbitals 
at the impurity sites and any of the six p-orbitals in the 
nearest-neighbor sites (notice that the hoppings between 
p levels between nearest neighbor and impurity sites are 
already contained in Eq.(31)). Finally, the Hund interac- 
tion between the localized spin and the spin of the hole 
in the (i-orbital at the impurity site is considered. If we 
were not working in the classical limit for the spins, a spin 
two operator (instead of S=5/2) will have to be used in 
Eqs.(34) and (31). 

Notice that the t sp hoppings are functions of the di- 
rection and, if the spin-orbit interaction is considered, 
of the quantum numbers j, rrij. Following the procedure 
described in Section II we know that in an fee lattice 



l(spa), E sy =m(spa), E sz =n(spcr). (35) 



Then the hoppings to the 12 neighbors are: 



-t£" = E sa {^v) = (-ir+ 



v (spcr) 
V2 



(36) 



if either fi = a or v = a, or zero otherwise. Notice that 
the minus (plus) sign in the hopping corresponds to the 
case in which \x and v have the same (opposite) sign. 

Applying the basis transformation given by Eqs.(12) 
and (13) to the p-orbitals the explicit form of the hop- 
pings between the d-orbital and the {j,mj) orbitals are 
obtained: 



= -t s 



(1+1) 

' V2 ' 



= -ts 



V2' 



P V2' dT '~3 P V6 ' 



" s p ns ' 



V6' 



,x,z _ . J_ ,y,z 



t x ' v — n t x - z — t y ' z — it 

c dT,5 ' dt,| ~~ dT,2 ~ sp V 3' 



,x,y _ ,x,z _ ,y,z _ ,x,y _ n 



f x,z =f y,z 



,x,z 

l Jt 1 



-it s 



1 .x,y _ _. (1+J) 

3' tsp V3 ' 



-t — t v ' z - t — 
sp V3' rfb-l SP V3' 



+x,y _ +x-z _ +y,z _ ,x,y _ n 

r di,l ~ r dl,§ - r dl,§ ~ l dl -\ - u ' 



,x,z _ ,y,z _ .. /2 x ,y _ _. (1 + 

W-k ~ ~ sp V 3' d U ~ sp V6 : 



di, 



<U>3 ~~ 



(i-0 



p V2 ' dl '-% ~ p v/2' 



"di,-* 



'"dj.-f _ S P /FT' ^ 



V2' 



V3 ' 



i 1 



Si _i — U ' S _I — S _1 ~ 4 ■ 



d| 



db-3 db 



(37) 



VI. NUMERICAL RESULTS IN FINITE 
SYSTEMS 

In this section, we address the results obtained from 
the study of Eq.(31) in a finite fee lattice. In recipro- 
cal space, the allowed momentum values form a cubic 
structure with 5ki = 2tt/ (otoiVj) where i — x, y, or z and 
Ni is the number of unit cubic cells in the sample along 
the i direction. Since in laboratory samples TV, is of the 
order of Avogadro's number we can replace fcj by a con- 
tinuous variable. However, this is not true in the small 
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FIG. 2: (color online) Band structure for GaAs obtained di- 
agonalizing Eq.(10) (continuous lines). The dashed lines in- 
dicate the results for the Luttinger-Kohn model. 



finite systems that can be studied numerically even with 
the most powerful computers currently available. Only 
when J = we can diagonalize Eq.(10) using continuous 
values of k. The results for GaAs along high symmetry 
directions in the first Brillouin zone (FBZ) are indicated 
by the continuous line in Fig. 2. As expected, the bottom 
of the band is at V and we observe the heavy hole and 
light hole bands along the high symmetry lines in the 
Brillouin zone shown in the figures. Since we are using 
"hole" language the top of the electronic valence band 
appears here as the bottom of the "hole" valence band. 
The dashed lines are the eigenvalues of Eq. ( A8) in Ref |7|, 
i.e., the Luttinger-Kohn model results. In order to check 
the agreement with our results around the point V we 
shifted our curves by the value 8tl x + At xx so that the 
bottom of our valence band is at zero. It can be seen that 
the agreement between the curves obtained with our tight 
binding model and the Luttinger-Kohn ones is excellent 
at the r point of the valence band. However, our model 
captures better the band dispersion away from the cen- 
ter of the Brillouin zone and the agreement with ARPES 
results is remarkable 4£ 

When J =/= and a finite number of magnetic impu- 
rities are considered, the diagonalization of Eq.(31) has 
to be performed in real space. In Fig. 3 we show the dis- 
crete energy eigenvalues along the high symmetry lines in 
the FBZ obtained numerically diagonalizing Eq.(31) us- 
ing finite clusters (still with J = 0). N ~ N x = N y = N z 
indicates the number of unit cells considered along each 
spatial direction. Numerically, we study clusters that 
contain N cubes of side ao along each of the 3 spa- 
tial directions. Since there are 4 ions associated to 
each site of the cube in an fee lattice, the total num- 
ber of Ga ions in the numerical simulation is given by 
Nca = iN 3 . This means that there is an equal num- 
ber of points inside the FBZ in momentum space. As 
already mentioned, the discrete lattice in momentum 



FIG. 3: (color online) Band structure for GaAs obtained di- 
agonalizing Eq.(10) (black lines). The symbols indicate re- 
sults obtained using finite lattices in real space diagonalizing 
Eq.(31) with J=0. As the lattice size increases, more momen- 
tum states in the First Brillouin zone become available. The 
number of momentum states available inside the FBZ is given 
by 4iV 3 . 



space is cubic. The side of the smallest cube is given by 
b = 2i:/Nao which is the size of the mesh corresponding 
to N cubic cells along each of the three spatial direc- 
tions. The first Brillouin zone has the shape of a trun- 
cated octahedron defined by well known high symmetry 
points such as L — (n/ao, 7r/ao, n/ao), X — (2ir/ao, 0, 0), 
U = (27r/a ,7r/2ao,7r/2a ), K = (37r/2a , 3tt/2o , 0), 



etc 



40 



The finiteness of the lattices imposes a constraint on 
the minimum Mn concentrations that can be reached: 
values of x larger than 0.2% can be studied in lattices 
with N = 4, i.e. 256 sites. The results presented in this 
work have been obtained using standard Monte Carlo 
simulations for the localized spins configurations which 
requires, at each iteration, the exact diagonalization of 
a AN 3 N x 4N 3 N Q fermionic matrix, where N Q = 6 for 
Eq.(31) and 8 when the "d-orbital" is included^ Using 
state-of-the-art computers such as the Cray XT3 super- 
computer at ORNL we are currently able to perform pro- 
duction runs for N — 4 and some runs with N = 5 and 6, 
and also study non cubic lattices, to monitor finite size 
effects. We are working on adapting newly developed 
(approximated) numerical techniques such as the TPEM 
approach^ to this problem with the hope of achieving 
production runs for = 6 and 8. 



A. Four Band Approximation 

In the case of III-V semiconductors with strong spin- 
orbit interaction, i.e., large A$o such as GaAs and GaSb, 
only the four j — 3/2 orbitals are often considered to 
study the top of the valence bandi&i This approximation 
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would certainly simplify numerical simulations since the 
number of degrees of freedom per lattice site is reduced 
from six to four. 

However, we have found an important drawback. As it 
can be seen in Fig. 4 the band dispersions under the four- 
orbital approximation (continuous line) agrees with the 
dispersions for the six-orbital case (dashed line) only in 
a very small region about T with radius 8ka^ w 0.02-7T 
while the study of many representative dopings such 
as x = 8.5% with p = 0.7 involves states at a radius 
Skao « 0.57T. In this situation the four-orbital approxi- 
mation no longer holds. Along the T — L, T — X, and 
r — U directions the heavy-hole dispersion is well repro- 
duced for a large momentum range but the very close in 
energy light-hole state is not captured by the 4 orbital ap- 
proximation. Even worse, along the T — K direction, the 
flatness characteristic of the heavy-hole band is missed in 
the 4 orbital approach. Thus, the 4 orbital model would 
only be adequate to describe very lightly doped cases 
(probably in the insulating regime) beyond the reach of 
finite lattice simulations and not very relevant for the 
high Tc regime. 




FIG. 4: (color online) Band structure for GaAs obtained di- 
agonalizing Eq.(10) (dashed lines). The band structure using 
the 4 orbital approximation is indicated with continuous lines. 



An example of the shortcomings of the 4 orbital model 
is the incorrect form of the magnetization M as a function 
of temperature (T) that it provides. In Fig. 5a the magne- 
tization curve for x = 8.5% Mn doped GaAs calculated 
with the 4-orbital approximation (circles) is compared 
with the one given by Eq.(31) (diamonds). The Curie- 
Weiss (CW) shape observed in experiments^ is only cap- 
tured by the 6-orbital model^l while the 4-orbital approx- 
imation reproduces the concave shape observed in previ- 
ous numerical efforts dealing mostly with single orbital 
modelsi^di In addition, the value of Tq is underestimated 
in this approach. 



B. Spin-Orbit Interaction 

As pointed out in the previous subsection, a prob- 
lem of early unbiased numerical efforts was the failure 
to reproduce the Curie- Weiss shape of the magnetiza- 
tion curves obtained in experiments of Mn doped GaAs 
with low compensation. While this shape was charac- 
teristic of most mean-field based approaches^, a con- 
cavity consistent with percolative behavior was observed 
otherwise 2*LL 

According to our results the CW behavior is very sen- 
sitive to the dispersion curves and the number of orbitals. 
Numerical simulations of the model presented here have 
shown that the spin-orbit interaction is crucial to obtain 
the CW shape of the magnetization curves.— In Fig. 5a 
we present magnetization curves obtained in lattices with 
256 Ga ions for J = 1.2eV, x = 8.5%, and the hole 
density p — 0.7. The curve with Aso — 0.34eV, indi- 
cated with diamonds, is in excellent agreement with the 
experimental data for Mn doped GaAs exhibiting CW 
behavior. It also provides an appropriate value of Tc- 
On the other hand, the curve with Aso = (indicated 
with squares), shows a linear increase of the magneti- 
zation below Tc- We have also used Eq.(31) to calcu- 
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FIG. 5: (color online) Magnetization vs temperature calcu- 
lated numerically with x=8.5% of Mn, and a density of holes 
p=0.72 for (a) GaAs; the diamonds indicate the results with 
spin-orbit interaction while the squares show the points for 
which the spin-orbit interaction has been neglected. Results 
using the 4-orbital approximation, i.e. Aso — » °o are de- 
noted by the circles, (b) GaSb (circles), GaAs (diamonds), 
GaP (squares), and GaN (triangles). The magnetization is 
measured as M — vM ■ M, with M the vectorial magnetiza- 
tion. As a consequence, for fully disordered spins, A4 is still 
nonzero due to the Si 2 = l contributions, causing a finite value 
at large temperatures (M(T — » oo) = 1/y/xNaa) unrelated 
to ferromagnetism. Thus, we plotted M = (A4 — M(T — » 
oo))/(l — A4(T — ► oo)), i.e. the background was subtracted. 



late the magnetization of other doped III-V compounds 
with higher and lower spin-orbit coupling than GaAs. 
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GaSb has a stronger spin-orbit coupling A50 = 0.75eV. 
Eq.(31) can be used for GaSb with the hoppings obtained 
using Eq.(7) and (15) using the Luttinger parameters 
(71,72,73) = (13.3, 4.4, 5. 7). 5 J is obtained using Eq.(30) 
and, since a = 6.101 5 , J=0.96eV. 

The magnetization, denoted by the circles in Fig. 5b 
shows CW behavior but Tc ~ 100K is lower than for 
GaAs at the same doping and compensation. The re- 
sult roughly agrees with mean-field estimates^ We also 
present results for Mn doped GaP which has a much 
smaller spin-orbit coupling (Ago = 0.08eV) than GaAs. 
The Luttinger parameters are (4.05, 0.49, 1.25)^ and we 
use J = 1.34eV. In this case we obtain a magnetiza- 
tion curve with an almost linear temperature dependence 
(squares in Fig. 5b). This linear behavior has been ob- 
tained in measurements of (Ga,Mn)P with 1.8% < x < 
3.8%£ 2 - but it is not clear whether it will continue up to 
8.5% doping if this value were reached. 

For completeness, we also present the results for 
(Ga,Mn)N. A Tc ~ 500K, well above room tempera- 
ture, is obtained in this case. This result is in agree- 
ment with mean-field predictions, assuming the doped 
Mn in the state Mn 2+ ,— i 53 i 54 and with some experimen- 
tal measurements 1 35 ' 46 The effects of Mn-ci-level partici- 
pation will be discussed in subsection E. 



C. J-Induced Splitting of the Valence Band at the 
T-Point 

Many studies of DMS are based on reasonable as- 
sumptions on the properties of the ground state. Some 
researchers work in the "high doping" regime assum- 
ing that the doped holes are uniformly distributed in 
the system and effectively doped into the valence band 
of the material. This is the so-called valence-band 
scenario^^ ' 56 i 57 in which the disorder configurations of 
the Mn ion are disregarded. The other approach, known 
as the impurity-band regime, 58 considers the limit of very 
low doping in which the holes electrically bound to the 
impurity cores start forming an impurity band due to the 
overlap of the wave functions. In this case the disordered 
positions of the Mn ions would play an important role. 
Although each scenario must be valid in opposite doping 
regimes, it is not clear when the crossover between the 
two occurs (since it is a function of x, J, and other param- 
eters) and which one better describes the experiments is 
controversial. Even for (Ga,Mn)As both descriptions are 
applied to the case of x=8.5%. 

A way of shedding light on this issue is by studying 
the evolution as a function of J, and for different val- 
ues of a;, of the four degenerate states at the top of the 
valence band (r point) of the parent compound. This 
splitting is due to the interaction of the doped holes with 
the localized spins and it was briefly described when the 
classical limit was discussed in section IV. A. For small 
values of J, at values of x in the metallic regime accord- 
ing to Mott's criterion, it is expected that the mean-field 



approach (assuming uniform hole distribution) should be 
a good approximation. The calculations indicate ^ 29 ! 30 
that the levels will split following a linear behavior in Jx, 
i.e., E ± Bq and E ± 3B G with B G given in Eq.(23) 
and (24). 

We have evaluated the eigenvalue distribution for dif- 
ferent Mn disorder configurations (indicated with differ- 
ent symbols in Fig. 6) as a function of J for x=8.5% 
(Fig. 6a) and 5% (Fig. 6b). For both values of x it can 
be seen that the mean-field prediction is accurate for 
J < 0.5eV; for J=1.2eV, i.e., the accepted value for 
GaAs, the eigenvalues have lower energies than the mean- 
field ones indicating a stronger distortion of the bands 
due to the interactions. However, the eigenvalues appear 
to be independent of the disorder configurations indicat- 
ing that the holes are non-localized; this uniform hole 
distribution continues up to J « 4eV which is the value 
at which our numerical simulations indicate that an im- 
purity band starts to develop in the density of states^ 
For larger values of J there is a large spread in the eigen- 
value distribution as a function of the disorder. Notice 
that the effects of disorder are much stronger for x—5% 
than for x—8.5% which corroborates the assumption that 
the role of disorder diminishes as x increases. It is inter- 
esting to notice that, in disagreement with the mean-field 
expectation that the spread in the eigenvalue distribution 
should increase with x, we actually observe that as x in- 
creases the width of the eigenvalues band decreases. 

■ 1 1 1 1 1 1 1 — 1 — 1 J...- 1 -. - H f ' I ' I ' I ' I ' I ' 1 - 



-9.5 - 
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FIG. 6: (color online) Numerically determined evolution of 
the four degenerate states at the top of the valence as a func- 
tion of J for (a) x = 0.085 and (b) for x — 0.05 in a 256 sites 
lattice. The different symbols indicate different Mn disorder 
configurations. The dotted lines correspond to the mean-field 
predictions. The insets show the range in J for which agree- 
ment between the numerical and mean-field results occur. 



These results indicate that J = 1.2eV corresponds to 
an intermediate region in which in the metallic regime 
the holes are not localized (if in the metallic regime of 
x)but the mean-field approximation needs corrections. 



13 



D. Dependence of Tc with x 



E. Coulomb Attraction at Mn Sites 



An important issue in DMS is the dependence of the 
Curie temperature with the amount of magnetic impuri- 
ties. Mean-field calculations 7 predict that Tc increases 
linearly with x but numerical simulations in simple sys- 
tems indicate that Tc reaches a maximum at an optimal 
value of x (x op t) and decreases afterwards.— The opti- 
mal value of x was found to be a function of J. For 
large J, x op t — > 1 while x op t ~ 0.25 for the values of 
J that capture the physics of (Ga,Mn)As. Experimen- 
tally, in (Ga,Mn)As, Tc has been found to increase with 
x up to about x=10%4 Recently, x=12.2%- 21.3% has 
been achieved^ but the highest Tc observed has been 
170K, i.e., no higher than the record value (T<7=173K) 
obtained with x=9%.— In fact, for both annealed and as 
grown samples it was observed that Tc reaches a broad 
maximum for x ~10%-15% and decreases afterwards. It 
is not clear whether the decrease is due to intrinsic prop- 
erties of the material or to a larger amount of compen- 
sation. In Fig. 7 we present results for Tc vs x obtained 
with Eq.(31) using the parameters for (Ga,Mn)As. Ex- 
perimental points are presented for comparison. As x 
increases longer thermalizations are needed in the nu- 
merical simulations. For p — 0.5, 0.75 and 1 we have 
obtained reasonable agreement with the experimental 
data for the annealed samples. We obtained a maximum 
T c = (225±20)K for x = 12.5% in the non-compensated 
case, i.e. p = 1. Thus, according to our calculations, 
even larger Mn doping of GaAs would not allow to reach 
a Curie temperature above room temperature. 
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FIG. 7: (color online) Numerically determined Tc as a func- 
tion of x for Ga i -jMiijAs for different values of compensation 
p. Experimental data are also shown. 



Numerical calculations of the density of states (DOS) 
using Eq.(31) for the case of (Ga,Mn)As have indicated 
that the chemical potential is in the valence band for 
dopings x=3% or higher^ 7 . As discussed in Section VA 
it has been argued by some researchers that this outcome 
may be the result of neglecting the effects of Coulomb at- 
traction between the localized Mn 2+ ions and the doped 
holes. In this subsection we will study the effects of 
adding the terms Eq.(32) and Eq.(33) to Eq.(31). In 
Ref|3J it was estimated that V — — 1.9eV for Mn doped 
GaAs. In Fig. 8 it can be seen that the DOS for both 
x=8.5% (Fig. 8a) and 3% (Fig. 8b) remains basically un- 
changed when either an onsite or a nearest neighbor 
Coulomb potential of the above magnitude is added. The 
chemical potential is still in the valence band and no hint 
of impurity band appears. Since even for the zero range 
potential the holes are already delocalized, there is no 
difference when the longer range potential is considered. 
Similar results (no IB in the DOS) are obtained for the 
case of (Ga,Mn)P using V=-2.4 eV. 34 This indicates that 
for the values of the parameters for (Ga,Mn)As, the ad- 
dition of Coulomb attraction is irrelevant in the range of 
doping of interest, i.e., x> 3%. 
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FIG. 8: (color online) Numerically determined density of 
states for Mn doped GaAs with (a) x — 0.085 and and 
p = 0.72 and (b) for x = 0.05 and p = 0.5. The different 
lines indicate different ranges of the Coulomb potential. Vnn 
indicates the value of the Coulomb potential at the impurity's 
nearest neighbor sites. 



However, for larger (but unphysical) values of J our 
results confirm the conclusions of Refl34l regarding the 
importance of considering a nearest-neighbor range po- 
tential in order to accurately obtain the value of x for 
which the crossover between the IB and VB description 
occurs. In Fig. 9 we show the DOS for the unphysical case 
in which the hopping parameters correspond to GaAs 
but we set J=7eV instead of the realistic value 1.2eV for 
x=8.5% (dashed line in Fig. 9). We chose J=7 eV because 
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it is close to the value for which a clearly separated IB 
develops in the DOSi^l By adding an on-site Coulomb 
potential of strength -3.5eV an IB develops (continuous 
line) and the chemical potential is located there. Clearly 
the holes are trapped at the Mn sites due to the large 
potential. However, if a nearest-neighbor range poten- 
tial with the same strength is applied, the holes become 
delocalized and the IB band in the DOS disappears (dot- 
dashed line). 
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FIG. 9: (color online) Numerically determined density of 
states for J = 7eV with x=8.5%. For V = (dashed line), on 
site V==-3.5 eV (continuous line), and nearest-neighbor range 
V==-3.5 eV (dot-dashed line). 



The above result indicates that studies of the crossover 
from impurity to valence band regimes performed with 
on-site rather than extended attractive potentials will 
overestimate the range of the impurity band regime. 



F. Mn d-orbitals 

Finally, in this subsection, we will consider the possi- 
bility of Mn 3+ doping in some III-V compounds such as 
GaN. 

In Fig. 10 the stars indicate the magnetization as a 
function of temperature obtained using the six-orbital 
Hamiltonian in Eq.(31) with the parameters for GaN 
on a 256 sites lattice which has already been presented 
in Fig. 5. For this material some experiment a 36 ' 38 and 
calculation o 37 i 61 ' 62 suggest that Mn d-levels may not be 
deep into the valence band. Thus, Mn 3+ , instead of 
Mn 2+ , may be present. To explore the consequences of 
this possibility we will use the 8-orbital model that results 
from the addition of Eq.(34) to Eq.(31). As mentioned in 
Section V.B the relative position of the d level Vd and the 
d — d and d — p hoppings will be input parameters. First, 
we will consider the case in which the hoppings involv- 
ing the d-orbitals are smaller than the hoppings among 
p-orbitals. Fixing tdd — — O.leV and t p d = O.leV, we 
present the magnetization vs T for various values of Vd 
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FIG. 10: (color online) Numerically determined magnetiza- 
tion vs T for Mn doped GaN with phenomenological d or- 
bitals. The open symbols correspond to tdd = —O.leV and 
tpd = O.leV' while the close symbols indicate tdd = 0.5eV and 
t v d = 0.5eV. The stars are data for the 6-orbital model. 



indicated in Fig. 10 (open symbols). For Vd — — 3.8eV, 
which corresponds to the (i-orbitals deep into the valence 
band (open circles) , the results are similar to the ones ob- 
tained for the 6-orbital case as expected. However, as Vd 
becomes more negative indicating that the <i-level moves 
into the gap it can be clearly seen that both the mag- 
netization and Tc decrease. However, it is important to 
notice that, even in the case of d-orbitals in the gap, the 
magnetization and Tc increase if the hoppings between 
d and p orbitals are allowed to be larger. For example, 
for | tpd | = | tdd | = 0.5eV and Vd = — 5.2eV, indicated 
with filled triangles in Fig. 10, the Curie temperature be- 
comes very close to room temperature. This behavior is 
in agreement with the ab initio results of Refl37l. 

Tc is related to the band structure as it can be seen in 
Fig. 11. As \Vd\ becomes larger an impurity band devel- 
ops in the DOS. When t p d << t pp the IB has very small 
dispersion indicating that the trapped holes are very lo- 
calized. When t p d increases the IB acquires dispersion 
and Tc becomes higher. 



VII. CONCLUSIONS 

We have presented a Hamiltonian that enables the nu- 
merical study of the physics of DMS with an equally ac- 
curate treatment of both the kinetic and magnetic inter- 
actions. Provided that the parent compound is a Zinc- 
Blende type semiconductor for which the Luttinger pa- 
rameters are known, there are no arbitrary parameters 
in the Hamiltonian. Only 6 orbitals per site have to be 
considered to obtain accurate Curie temperatures and 
magnetization curves of materials where Mn 2+ is present, 
which makes possible to perform numerical simulations 
in large systems with present day computers. Besides re- 
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FIG. 11: (color online) Numerically determined DOS for Mn 
doped GaN with phenomenological d orbitals. The continu- 
ous lines correspond to tdd = — O.leV and t p d = O.leV for 
the values of Vd indicated while the dashed line indicates 
tdd = 0.5eV and t p d = .5eV. Data for the 6-orbital model 
are denoted with dashed-dotted line. 



producing well established experimental results for GaAs 
with less than 10% Mn doping, the higher doping regime 
was explored. While naively it could have been expected 
that Tq would increase continuously with x, we observed 
that it reaches a maximum value at an intermediate value 
of the doping. For GaAs we found that the highest value 
of Tc, approximately 220K, occurs x—12.5% and no com- 
pensation, indicating that higher doping of GaAs will not 
lead to Tp's above room temperature. The doping for 
which the maximum Tq occurs appears to be in general 
agreement with recent experimental results. 

The possibility of Mn 3+ doping in GaP or GaN has 
been considered phenomenologically. We found that Tc 
will be very reduced if the d levels are located at energies 
much higher than the top of the valence band and if the 
pd hybridization is weak. However, if the hybridization 
is strong, Mn doped GaN could be a good candidate for 
high T c . 

It has also been shown that the Coulomb attraction be- 
tween the hole and the Mn ions expected at low doping 
does not play a role for x larger than 3% in GaAs. As a 
consequence, we have not observed an impurity band that 
should be present at lower dopings or even in the metallic 
regime if the magnetic interaction were stronger. The rel- 
ative weakness of the magnetic interaction was confirmed 
by studying the splitting of the top of the valence band as 
a function of J. For values of J within the experimental 
range for Mn doped GaAs, we have observed a departure 



from the mean-field predictions. However, the holes do 
not appear to be localized at the Mn impurity sites since 
the eigenvalue distribution is independent of the particu- 
lar disorder configuration. Localization is observed only 
at unrealistically high values of J. 

Summarizing, we have shown that the real space 
Hamiltonian presented here reproduces the top of the 
valence band given by the Luttinger-Kohn model for the 
undoped parent compounds and, as such, it is guaran- 
teed to provide an excellent framework for the numerical 
study of lightly doped Zinc-blende type magnetic semi- 
conductors. 
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IX. APPENDIX I 

Here we provide the change of basis matrices used in 
Eq.(12) and (13). 
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